close all


% data summary
%(10, 9): 
% {'20230804_100VR_5G_10kHz_K2_scan_2',  '20230808_100VR_5G_10kHz_K2_scan_1'}; 
% cyclecut =  [130 149; 530 546];

%(10, 10): 
% {'20230812_100VR_5G_10kHz_K2_scan_1','20230807_100VR_5G_10kHz_K2_scan_2'};
% cyclecut = [];

%(9, 10): 
% {'20230805_100VR_5G_10kHz_K2_scan_1','20230806_100VR_5G_10kHz_K2_scan_1','20230803_100VR_5G_10kHz_K2_scan_8'}; 
% cyclecut = [396 722];

%(9, 9):
% {'20230810_100VR_5G_10kHz_K2_scan_2','20230809_100VR_5G_10kHz_K2_scan_1'}; 
% cyclecut = [];

% background (10,10) at 50G
% '20230814_100VR_50G_10kHz_K2_scan_1'

Files = {{'20230804_100VR_5G_10kHz_K2_scan_2',  '20230808_100VR_5G_10kHz_K2_scan_1'} ...
     {'20230812_100VR_5G_10kHz_K2_scan_1','20230807_100VR_5G_10kHz_K2_scan_2'}, ...
     {'20230810_100VR_5G_10kHz_K2_scan_2','20230809_100VR_5G_10kHz_K2_scan_1'} ...
     {'20230805_100VR_5G_10kHz_K2_scan_1','20230806_100VR_5G_10kHz_K2_scan_1','20230803_100VR_5G_10kHz_K2_scan_8'}, ...    
     {'20230814_100VR_50G_10kHz_K2_scan_1'}};
Cycle = {[130 149; 530 546],[],[],[396 722], []};
Names = {[10,9],[10,10],[9,9],[9,10], [10,10]};

%%
%Input parameters

h4 = figure(); defFontSize = 16;
figure(h4);
% set(h4, 'Position', [scrsz(3) 1 scrsz(3) scrsz(4)])


for ss = 1:length(Files)

filename = Files{ss};{'20230814_100VR_50G_10kHz_K2_scan_1'};
DeltaT = 45*1e3; %units are ns. Represents time between ODT off and REMPI/UV pulse
%DeltaTlist = DeltaT;
DeltaTlist = linspace(0,DeltaT,80);
shotNumStart = 000; % 300
maxShot = 2*10000;

%Changes for decoherance
shotNumEnd = maxShot./2;
Momemtum_filter_scale = 1;


nXmomfilter = 3; %determines how many sigmas to use for filtering criterion. 1 is 1-sigma
nYmomfilter = 3; %determines how many sigmas to use for filtering criterion. 1 is 1-sigma
nTOFfilter = 3; %determines how many sigmas to use for filtering criterion. 1 is 1-sigma
Bfield = 5; %Value of the Bfield in Gauss

state_pair =  Names{ss};
cyclesCut = Cycle{ss};

%%
%Flags for different functions
TOFfilterFlag = 1;
XmomfilterFlag = 1;
YmomfilterFlag = 1;
VariableScanFlag = 0;

CoinShotDistribution = 0;
shot_bin_size = 10000;

%%
TOFAllShots = [];
xPosAllShots = [];
yPosAllShots = [];
shotNumAllShots = [];
xVarAllShots1 = [];
xVarAllShots2 = [];
flagREMPILockAllShots = [];
for i = 1:length(filename)
    filepath = [folderpath, filename{i}];
    if (~exist([filepath,'.dat'],'file'))
        error([filepath,'.dat',' needs to be created using mainMCPtxtAnalysisTiming.m first!'])
    else
        realHit = dlmread([filepath,'.dat']);
    end
    TOFAllShotsSub = realHit(:,1);
    xPosAllShotsSub = realHit(:,2);
    yPosAllShotsSub = realHit(:,3);
    shotNumAllShotsSub = realHit(:,4);
    xVarAllShots1Sub = realHit(:,5);
    xVarAllShots2Sub = realHit(:,6);
    flagREMPILockAllShotsSub = realHit(:,7);
    TOFAllShots = [TOFAllShots;TOFAllShotsSub];
    xPosAllShots = [xPosAllShots;xPosAllShotsSub];
    yPosAllShots = [yPosAllShots;yPosAllShotsSub];
    shotNumAllShots = [shotNumAllShots;shotNumAllShotsSub];
    xVarAllShots1 = [xVarAllShots1;xVarAllShots1Sub];
    xVarAllShots2 = [xVarAllShots2;xVarAllShots2Sub];
    flagREMPILockAllShots = [flagREMPILockAllShots;flagREMPILockAllShotsSub];
end
%%
%Define some relevant constants
mK2amu = 2*(39.96399848); %mass of K2 in amu
mRb2amu = 2*(86.909180520); %mass of Rb2 in amu
A = 15.23; %in mm/sqrt(cm^-1/V)
VR = 99; %volts


%%

    % Calculate the total number of experimental cycles
    totCycleNum = 1;
    cycleNumAllShots = zeros(length(realHit),1);
    cycleNumAllShots(1) = 1;
    
    for i = 2:length(shotNumAllShots)
        if shotNumAllShots(i) < shotNumAllShots(i-1)
            totCycleNum = totCycleNum + 1;
        end
        cycleNumAllShots(i) = totCycleNum;
    end
    
%     cyclesCut = [1100 1206];[600 1066];[300 2000]; [635 1210]; %[1 429;1160 1190]; % [100 200;300 400]
    [numCutPair,~] = size(cyclesCut);
    % Comments about "cyclesCut"
    % Each pair of values specify an interval of cycles to be cut out of the data analysis
    % Example: to cut from shots 100 - 200 and 300 - 400, use the format cyclesCut = [100 200;300 400];
    % To not cut anything, put [0 0]
    cycleNumInd = 1:length(cycleNumAllShots);
    for i = 1:numCutPair
        tmp = find((cycleNumAllShots <= cyclesCut(i,1))|(cycleNumAllShots >= cyclesCut(i,2))); % find indices of shots that will survive the cyclesCut
        cycleNumInd = intersect(cycleNumInd,tmp);
    end


    
    %% Define data vectors filtered by shot numbers and cycle numbers
    shotNumIndAll = find((shotNumAllShots >= shotNumStart)&(shotNumAllShots <= shotNumEnd));
    
    filterNumIndAll = intersect(shotNumIndAll, cycleNumInd);
    
    TOFAll = TOFAllShots(filterNumIndAll);
    xPosAll = xPosAllShots(filterNumIndAll);
    yPosAll = yPosAllShots(filterNumIndAll);
    shotNumAll = shotNumAllShots(filterNumIndAll);
    cycleNumAll = cycleNumAllShots(filterNumIndAll);
    xVarAll1 = xVarAllShots1(filterNumIndAll);
    xVarAll2 = xVarAllShots2(filterNumIndAll);
    flagREMPILockAll = flagREMPILockAllShots(filterNumIndAll);
    flagREMPILockAll = ones(length(flagREMPILockAll),1);
    
    LockedIndex = find(flagREMPILockAll == 1);
    LockedCycleNumAll = cycleNumAll(LockedIndex);
    LockedCycleNum = length(unique(LockedCycleNumAll));
    
    speciesNames = {'K','K_2','Rb','KRb','K_2Rb','Rb_2','KRb_2','K_2Rb_2'};
    
    speciesMasses = [40,80,87,127,167,174,214,254];
    
    TOFs = [48930, 69196, 72150, 87180, 99971, 102045, 113167, 122000]*1.096 - 19900; 
    TOFWindow = 1800;
    
    %% ROI window parameters for 99 V and 30 G
    theta = 45;
    if Bfield == 30
        xCentROISpecies = [12.95, 5.024, 5.05, 1.37, 0, -0.72, 0, 0]; %30 Gauss
        yCentROISpecies = [-0.85, -0.5887, -0.85, -0.85, 0, -0.85, 0, 0]; %30 Gauss
    elseif Bfield == 5
        xCentROISpecies = [-11.181, -12.2, -10.9494, -11.333, -11.55, -11.55, -11.55, -11.55]; % 5 Gauss
        yCentROISpecies = [-1.2188, -0.623, -0.5559, -0.8894, -0.8, -0.8, -0.75, -0.75]; % 5 Gauss
    end
    windowROISpecies = [10, 10, 10, 20, 10, 10, 10, 5];
    windowROISpecies = [40, 40, 40, 40, 40, 40, 40, 40]*2;
    xMinROISpecies = xCentROISpecies - windowROISpecies./2;
    xMaxROISpecies = xCentROISpecies + windowROISpecies./2;
    yMinROISpecies = yCentROISpecies - windowROISpecies./2;
    yMaxROISpecies = yCentROISpecies + windowROISpecies./2;
    
    xPosAllRot = cosd(theta)*xPosAll - sind(theta)*yPosAll;
    yPosAllRot = sind(theta)*xPosAll + cosd(theta)*yPosAll;
    
    K2Ind = find((TOFAll >= TOFs(2) - TOFWindow/2) & (TOFAll <= TOFs(2) + TOFWindow/2) ...
        & (xPosAllRot >= xMinROISpecies(2)) & (xPosAllRot <= xMaxROISpecies(2)) ...
        & (yPosAllRot >= yMinROISpecies(2)) & (yPosAllRot <= yMaxROISpecies(2))); 
%         & (flagREMPILockAll == 1));
    Rb2Ind = find((TOFAll >= TOFs(6) - TOFWindow/2) & (TOFAll <= TOFs(6) + TOFWindow/2) ...
        & (xPosAllRot >= xMinROISpecies(6)) & (xPosAllRot <= xMaxROISpecies(6)) ...
        & (yPosAllRot >= yMinROISpecies(6)) & (yPosAllRot <= yMaxROISpecies(6))); 
%         & (flagREMPILockAll == 1));
    
    
    TOFK2 = TOFAll(K2Ind);
    xPosK2 = xPosAll(K2Ind);
    yPosK2 = yPosAll(K2Ind);
    shotNumK2 = shotNumAll(K2Ind);
    cycleNumK2 = cycleNumAll(K2Ind);
    meanTOFK2 = mean(TOFK2);
    meanTOFK2Offset = meanTOFK2 + 29944.7;
    
    TOFRb2 = TOFAll(Rb2Ind);
    xPosRb2 = xPosAll(Rb2Ind);
    yPosRb2 = yPosAll(Rb2Ind);
    shotNumRb2 = shotNumAll(Rb2Ind);
    cycleNumRb2 = cycleNumAll(Rb2Ind);
    meanTOFRb2 = mean(TOFRb2);
    meanTOFRb2Offset = meanTOFRb2 + 29944.7;
    
    COMFitFunction = fittype('a1*(1-exp(-x/a9))+a2','independent','x',...
        'coefficients',{'a1','a2','a9'});
        
    % COMFitxPosK2 = fit(shotNumK2,xPosK2,COMFitFunction,'Start',[-2.638 5.601 171.6]);

    x_mean_K2 = median(xPosK2);
    y_mean_K2 = median(yPosK2);
    x_mean_Rb2 = median(xPosRb2);
    y_mean_Rb2 = median(yPosRb2);
    filter_range =2;
    ind_to_fit_K2_x = find(xPosK2>=x_mean_K2-filter_range & xPosK2<=x_mean_K2+filter_range);
    ind_to_fit_K2_y = find(yPosK2>=y_mean_K2-filter_range & yPosK2<=y_mean_K2+filter_range);
    ind_to_fit_Rb2_x = find(xPosRb2>=x_mean_Rb2-filter_range & xPosRb2<=x_mean_Rb2+filter_range);
    ind_to_fit_Rb2_y = find(yPosRb2>=y_mean_Rb2-filter_range & yPosRb2<=y_mean_Rb2+filter_range);



 
    % COMFitxPosK2 = fit(shotNumK2,xPosK2,COMFitFunction,'Start',[-2.638 20 171.6]);
      COMFitxPosK2 = fit(shotNumK2(ind_to_fit_K2_x),xPosK2(ind_to_fit_K2_x),COMFitFunction,'Start',[-2.638 5.601 184]);
    pointsEvalCOMFitxPosK2 = linspace(shotNumStart,shotNumEnd,shotNumEnd-shotNumStart+1);
    EvalCOMFitxPosK2 = COMFitxPosK2(pointsEvalCOMFitxPosK2);
    % COMFityPosK2 = fit(shotNumK2,yPosK2,COMFitFunction,'Start',[4.688 -8.698 184.7]);
    % COMFityPosK2 = fit(shotNumK2,yPosK2,COMFitFunction,'Start',[4.688 2 184.7]);
    COMFityPosK2 = fit(shotNumK2(ind_to_fit_K2_y),yPosK2(ind_to_fit_K2_y),COMFitFunction,'Start',[4.688 2 184.7]);
    pointsEvalCOMFityPosK2 = linspace(shotNumStart,shotNumEnd,shotNumEnd-shotNumStart+1);
    EvalCOMFityPosK2 = COMFityPosK2(pointsEvalCOMFityPosK2);
    RelxPosK2 = xPosK2 - COMFitxPosK2(shotNumK2);
    RelyPosK2 = yPosK2 - COMFityPosK2(shotNumK2);
    
    % COMFitxPosRb2 = fit(shotNumRb2,xPosRb2,COMFitFunction,'Start',[-2.177 1.106 143.9]);
    COMFitxPosRb2 = fit(shotNumRb2(ind_to_fit_Rb2_x),xPosRb2(ind_to_fit_Rb2_x),COMFitFunction,'Start',[3 5 184.7]);
    % COMFitxPosRb2 = fit(shotNumRb2,xPosRb2,COMFitFunction,'Start',[-7.8 0 143.9]);
%     COMFitxPosRb2 = fit(shotNumRb2,xPosRb2,COMFitFunction,'Start',[-7.8 15 143.9]);
    pointsEvalCOMFitxPosRb2 = linspace(shotNumStart,shotNumEnd,shotNumEnd-shotNumStart+1);
    EvalCOMFitxPosRb2 = COMFitxPosRb2(pointsEvalCOMFitxPosRb2);
    % COMFityPosRb2 = fit(shotNumRb2,yPosRb2,COMFitFunction,'Start',[3.407 -3.408 170.7]);
    % COMFityPosRb2 = fit(shotNumRb2,yPosRb2,COMFitFunction,'Start',[3.407 1 170.7]);
    COMFityPosRb2 = fit(shotNumRb2(ind_to_fit_Rb2_y),yPosRb2(ind_to_fit_Rb2_y),COMFitFunction,'Start',[3.407 1 184]);
    pointsEvalCOMFityPosRb2 = linspace(shotNumStart,shotNumEnd,shotNumEnd-shotNumStart+1);
    EvalCOMFityPosRb2 = COMFityPosRb2(pointsEvalCOMFityPosRb2);
    RelxPosRb2 = xPosRb2 - COMFitxPosRb2(shotNumRb2);
    RelyPosRb2 = yPosRb2 - COMFityPosRb2(shotNumRb2);
    
    %%
    shotCycleNumK2 = [cycleNumK2,shotNumK2];
    shotCycleNumRb2 = [cycleNumRb2,shotNumRb2];
    [shotCycleCoin,K2CoinInd,Rb2CoinInd] = intersect(shotCycleNumK2,shotCycleNumRb2,'rows');
    prelimCoinTOFK2 = TOFK2(K2CoinInd);
    prelimCoinshotNumK2 = shotNumK2(K2CoinInd);
    prelimCoincycleNumK2 = cycleNumK2(K2CoinInd);
    prelimCoinRelxPosK2 = RelxPosK2(K2CoinInd);
    prelimCoinRelyPosK2 = RelyPosK2(K2CoinInd);
    prelimCoinTOFRb2 = TOFRb2(Rb2CoinInd);
    prelimCoinshotNumRb2 = shotNumRb2(Rb2CoinInd);
    prelimCoincycleNumRb2 = cycleNumRb2(Rb2CoinInd);
    prelimCoinRelxPosRb2 = RelxPosRb2(Rb2CoinInd);
    prelimCoinRelyPosRb2 = RelyPosRb2(Rb2CoinInd);
    prelimCoinCounts = length(prelimCoinRelxPosK2);
    
    NumShiftPts = 50;
    ShiftVector = -NumShiftPts:1:NumShiftPts;
    EqualShotIndex = ((length(ShiftVector) - 1)/2) + 1;
    ShiftVector = ShiftVector';
    CoinCountsShiftNoFilter = zeros(length(ShiftVector),1);
    CoinCountsShiftMomFilter = zeros(length(ShiftVector),1);
    CoinCountsShiftTOFFilter = zeros(length(ShiftVector),1);
    CoinCountsShiftAllFilter = zeros(length(ShiftVector),1);
    sigmaXmomfilter = 0.41*0.56/Momemtum_filter_scale; %resolution in mm, was 0.41*0.56
    sigmaYmomfilter = 0.41*0.56/Momemtum_filter_scale; %resolution in mm, was 0.41*0.56
    sigmaTOFfilter = 20.5*0.56; %resolution in ns
    L1 = 0.041;   %[m] length of VMI optics
    L2 = 1 - L1;  %[m] E field free TOF tube length
    eta = ((L2/(2*L1))^2-1);
    RatioK = (eta/meanTOFK2Offset)*DeltaTlist - 1;
    RatioRb = (eta/meanTOFRb2Offset)*DeltaTlist - 1;
    RatioKSquared = RatioK.*RatioK;
    RatioRbSquared = RatioRb.*RatioRb;
    TOFcomparison = nTOFfilter*sigmaTOFfilter*sqrt(RatioKSquared+RatioRbSquared);
    for i = 1:length(ShiftVector)
        shotNumK2Shift = shotNumK2 + ShiftVector(i);
        shotCycleNumK2Shift = [cycleNumK2,shotNumK2Shift];
        [shotCycleCoinShift,K2CoinIndShift,Rb2CoinIndShift] = intersect(shotCycleNumK2Shift,shotCycleNumRb2,'rows');
        prelimCoinTOFK2Shift = TOFK2(K2CoinIndShift);
        prelimCoinRelxPosK2Shift = RelxPosK2(K2CoinIndShift);
        prelimCoinRelyPosK2Shift = RelyPosK2(K2CoinIndShift);
        prelimCoinTOFRb2Shift = TOFRb2(Rb2CoinIndShift);
        prelimCoinRelxPosRb2Shift = RelxPosRb2(Rb2CoinIndShift);
        prelimCoinRelyPosRb2Shift = RelyPosRb2(Rb2CoinIndShift);
        CoinCountsShiftNoFilter(i) = length(prelimCoinRelxPosK2Shift);
        CoinIndXmomfilterShift = find(abs(prelimCoinRelxPosK2Shift + sqrt(mRb2amu/mK2amu)*prelimCoinRelxPosRb2Shift)...
            <= nXmomfilter*sigmaXmomfilter*sqrt(1+(mRb2amu/mK2amu)));
        CoinIndYmomfilterShift = find(abs(prelimCoinRelyPosK2Shift + sqrt(mRb2amu/mK2amu)*prelimCoinRelyPosRb2Shift)...
            <= nYmomfilter*sigmaYmomfilter*sqrt(1+(mRb2amu/mK2amu)));
        MomentumFilterIndexShift = intersect(CoinIndXmomfilterShift,CoinIndYmomfilterShift);
        CoinCountsShiftMomFilter(i) = length(prelimCoinRelxPosK2Shift(MomentumFilterIndexShift));
        prelimCoinRelTOFK2Shift = prelimCoinTOFK2Shift - meanTOFK2;
        prelimCoinRelTOFRb2Shift = prelimCoinTOFRb2Shift - meanTOFRb2;
        WeightedprelimCoinRelTOFK2Shift = prelimCoinRelTOFK2Shift*RatioRb;
        WeightedprelimCoinRelTOFRb2Shift = prelimCoinRelTOFRb2Shift*RatioK;
        TotalprelimCoinRelTOFMatrixShift = abs(WeightedprelimCoinRelTOFRb2Shift + WeightedprelimCoinRelTOFK2Shift);
        TOFDifferenceShift = TotalprelimCoinRelTOFMatrixShift - TOFcomparison;
        TOFConditionalMatrixShift = TOFDifferenceShift <= 0;
        TOFConditionalVectorShift = sum(TOFConditionalMatrixShift,2);
        CoinIndTOFfilterShift = find(TOFConditionalVectorShift > 0);
        CoinCountsShiftTOFFilter(i) = length(prelimCoinRelxPosK2Shift(CoinIndTOFfilterShift));
        FullFilterIndexShift = intersect(CoinIndTOFfilterShift,MomentumFilterIndexShift);
        CoinCountsShiftAllFilter(i) = length(prelimCoinRelxPosK2Shift(FullFilterIndexShift));
        if i == EqualShotIndex
            CoinRelxPosK2 = prelimCoinRelxPosK2Shift(FullFilterIndexShift);
            CoinRelyPosK2 = prelimCoinRelyPosK2Shift(FullFilterIndexShift);
            CoinRelxPosRb2 = prelimCoinRelxPosRb2Shift(FullFilterIndexShift);
            CoinRelyPosRb2 = prelimCoinRelyPosRb2Shift(FullFilterIndexShift);
            prelimCoinCountsindexvector = 1:1:length(prelimCoinRelxPosK2Shift);
            NonCoinCountsIndex = setdiff(prelimCoinCountsindexvector,FullFilterIndexShift);
            NotCoinRelxPosK2 = prelimCoinRelxPosK2Shift(NonCoinCountsIndex);
            NotCoinRelyPosK2 = prelimCoinRelyPosK2Shift(NonCoinCountsIndex);
            NotCoinRelxPosRb2 = prelimCoinRelxPosRb2Shift(NonCoinCountsIndex);
            NotCoinRelyPosRb2 = prelimCoinRelyPosRb2Shift(NonCoinCountsIndex);
        end
    end
    
    BkdCoinCountsShiftNoFilter = CoinCountsShiftNoFilter;
    BkdCoinCountsShiftMomFilter = CoinCountsShiftMomFilter;
    BkdCoinCountsShiftTOFFilter = CoinCountsShiftTOFFilter;
    BkdCoinCountsShiftAllFilter = CoinCountsShiftAllFilter;
    
    BkdCoinCountsShiftNoFilter(EqualShotIndex) = [];
    BkdCoinCountsShiftMomFilter(EqualShotIndex) = [];
    BkdCoinCountsShiftTOFFilter(EqualShotIndex) = [];
    BkdCoinCountsShiftAllFilter(EqualShotIndex) = [];
    
    MeanBkdNoFilter = round(mean(BkdCoinCountsShiftNoFilter),1);
    MeanBkdMomFilter = round(mean(BkdCoinCountsShiftMomFilter),1);
    MeanBkdTOFFilter = round(mean(BkdCoinCountsShiftTOFFilter),1);
    MeanBkdAllFilter = round(mean(BkdCoinCountsShiftAllFilter),1);
    STDBkdNoFilter = round(std(BkdCoinCountsShiftNoFilter),1);
    STDBkdMomFilter = round(std(BkdCoinCountsShiftMomFilter),1);
    STDBkdTOFFilter = round(std(BkdCoinCountsShiftTOFFilter),1);
    STDBkdAllFilter = round(std(BkdCoinCountsShiftAllFilter),1);
    PeakErrorNoFilter = round(sqrt(CoinCountsShiftNoFilter(EqualShotIndex)),1);
    PeakErrorMomFilter = round(sqrt(CoinCountsShiftMomFilter(EqualShotIndex)),1);
    PeakErrorTOFFilter = round(sqrt(CoinCountsShiftTOFFilter(EqualShotIndex)),1);
    PeakErrorAllFilter = round(sqrt(CoinCountsShiftAllFilter(EqualShotIndex)),1);
    DifferenceErrorNoFilter = round(1000*sqrt(CoinCountsShiftNoFilter(EqualShotIndex) + STDBkdNoFilter^2)/LockedCycleNum,1);
    DifferenceErrorMomFilter = round(1000*sqrt(CoinCountsShiftMomFilter(EqualShotIndex) + STDBkdMomFilter^2)/LockedCycleNum,1);
    DifferenceErrorTOFFilter = round(1000*sqrt(CoinCountsShiftTOFFilter(EqualShotIndex) + STDBkdTOFFilter^2)/LockedCycleNum,1);
    DifferenceErrorAllFilter = round(1000*sqrt(CoinCountsShiftAllFilter(EqualShotIndex) + STDBkdAllFilter^2)/LockedCycleNum,1);
   

        
    %% Plot of coincidence counts after filtering

    
    subplot(3,2,ss)
    plot(ShiftVector,CoinCountsShiftAllFilter,'b-o','LineWidth',2,'MarkerSize',8)
    % plot(ShiftVector,CoinCountsShiftAllFilter/LockedCycleNum*1000,'b-o','LineWidth',2,'MarkerSize',8)
    % round(1000*(CoinCountsShiftAllFilter(EqualShotIndex)-MeanBkdAllFilter)/LockedCycleNum,1)
    hold on
    % plot(ShiftVector,MeanBkdAllFilter*ones(length(ShiftVector),1)/LockedCycleNum*1000,'r--','LineWidth',3)
     plot(ShiftVector,MeanBkdAllFilter*ones(length(ShiftVector),1),'r--','LineWidth',3)
    hold off
    xlim([min(ShiftVector), max(ShiftVector)])
    % ylim([0, 1.2*max(CoinCountsShiftAllFilter/LockedCycleNum*1000)])
    ylim([0, 1.2*max(CoinCountsShiftAllFilter)])
    xlabel('Shot number difference', 'FontSize', 12)
    ylabel('Ion counts per 1000 cycles', 'FontSize', 12)
    text(-10, 1.1*max(CoinCountsShiftAllFilter), {['{\it N}_{0} = ',num2str(CoinCountsShiftAllFilter(EqualShotIndex)),' ',char(177),' ',num2str(PeakErrorAllFilter/LockedCycleNum*1000)]}, 'FontSize', 12, 'Color','b')
    %text(0.4*max(ShiftVector), 1.15*max(CoinCountsShiftAllFilter), {['Total Cycles = ',num2str(LockedCycleNum)]}, 'FontSize', 18, 'Color','k')
    text(13, 1.*max(CoinCountsShiftAllFilter),...
        {['Total Cycles = ',num2str(LockedCycleNum)],'{\it N}_{coin} = {\it N}_{0} - {\it N}_{bkgd}',['          = ',num2str(round(1000*(CoinCountsShiftAllFilter(EqualShotIndex)-MeanBkdAllFilter)/LockedCycleNum,1)),' ',char(177),' ',num2str(DifferenceErrorAllFilter)],...
        '      (cts/1000 cycles)'},...
        'FontSize', 12, 'Color','k')
    text(min(ShiftVector)+2, 0.3*(1.*max(CoinCountsShiftAllFilter/LockedCycleNum*1000) - MeanBkdAllFilter), {['{\it N}_{bkgd} = ',num2str(MeanBkdAllFilter),' ',char(177),' ',num2str(STDBkdAllFilter/LockedCycleNum*1000)]}, 'FontSize', 12, 'Color','r')
    % text(min(ShiftVector)+2, 1.05*max(CoinCountsShiftAllFilter), {'TOF + Momentum Filtering',['n_{Mom} = ',num2str(nXmomfilter)],['\sigma_{Mom} = ',num2str(sigmaXmomfilter),' mm'],['n_{TOF} = ',num2str(nTOFfilter)],['\sigma_{TOF} = ',num2str(sigmaTOFfilter),' ns']}, 'FontSize', 12, 'Color','k')
    title(['State Pair = (',num2str(state_pair(1)),',',num2str(state_pair(2)),')'],'fontweight','bold','FontSize', 12)
    set(gca,'fontsize',12)
    
    set(gcf,'color','white')
end

